Stan

MKTG 585R Quantitative Marketing Pre-PhD Seminar

Probabilistic Programming

Stan is a probabilistic programming language that provides a general-purpose sampler using Hamiltonian Monte Carlo.

  • You specify your model in terms of probability functions.
  • Stan automatically compiles a C++ HMC sampler for your model.

Stan(islaw) Ulam

When we use ulam(), {rethinking} is translating our model into Stan, which then compiles an HMC sampler in C++. We want to start writing in Stan code directly, and the ulam() translation will help us get there.

Let’s go back to the foxes data.

# Estimate the direct causal effect of avgfood on weight.
fit <- ulam(
  alist(
    weight ~ dnorm(mu, sigma),
    mu <- beta0 + beta_food * avgfood + beta_group * groupsize,
    beta0 ~ dnorm(0, 0.2),
    c(beta_food, beta_group) ~ dnorm(0, 0.5),
    sigma ~ dexp(1)
  ), 
  data = foxes_list, # Specify the data list instead of a data frame.
  chains = 4,        # Specify the number of chains.
  cores = 4,         # Specify the number of cores to run in parallel.
  log_lik = TRUE,    # To compute model fit via WAIC and PSIS.
  cmdstan = TRUE     # Specify cmdstan = TRUE to use cmdstanr instead of rstan.
)
In file included from src/cmdstan/main.cpp:1:
In file included from src/cmdstan/command.hpp:4:
In file included from src/cmdstan/arguments/arg_data.hpp:4:
In file included from src/cmdstan/arguments/categorical_argument.hpp:4:
In file included from src/cmdstan/arguments/argument.hpp:4:
In file included from stan/src/stan/callbacks/writer.hpp:4:
In file included from stan/lib/stan_math/lib/boost_1.78.0/boost/lexical_cast.hpp:31:
In file included from stan/lib/stan_math/lib/boost_1.78.0/boost/lexical_cast/bad_lexical_cast.hpp:28:
In file included from stan/lib/stan_math/lib/boost_1.78.0/boost/throw_exception.hpp:24:
stan/lib/stan_math/lib/boost_1.78.0/boost/assert/source_location.hpp:75:14: warning: 'sprintf' is deprecated: This function is provided for compatibility reasons only.  Due to security concerns inherent in the design of sprintf(3), it is highly recommended that you use snprintf(3) instead. [-Wdeprecated-declarations]
        std::sprintf( buffer, ":%ld", static_cast<long>( line() ) );
             ^
/Applications/Xcode.app/Contents/Developer/Platforms/MacOSX.platform/Developer/SDKs/MacOSX.sdk/usr/include/stdio.h:188:1: note: 'sprintf' has been explicitly marked deprecated here
__deprecated_msg("This function is provided for compatibility reasons only.  Due to security concerns inherent in the design of sprintf(3), it is highly recommended that you use snprintf(3) instead.")
^
/Applications/Xcode.app/Contents/Developer/Platforms/MacOSX.platform/Developer/SDKs/MacOSX.sdk/usr/include/sys/cdefs.h:215:48: note: expanded from macro '__deprecated_msg'
        #define __deprecated_msg(_msg) __attribute__((__deprecated__(_msg)))
                                                      ^
In file included from src/cmdstan/main.cpp:1:
In file included from src/cmdstan/command.hpp:4:
In file included from src/cmdstan/arguments/arg_data.hpp:4:
In file included from src/cmdstan/arguments/categorical_argument.hpp:4:
In file included from src/cmdstan/arguments/argument.hpp:4:
In file included from stan/src/stan/callbacks/writer.hpp:4:
In file included from stan/lib/stan_math/lib/boost_1.78.0/boost/lexical_cast.hpp:31:
In file included from stan/lib/stan_math/lib/boost_1.78.0/boost/lexical_cast/bad_lexical_cast.hpp:28:
In file included from stan/lib/stan_math/lib/boost_1.78.0/boost/throw_exception.hpp:24:
stan/lib/stan_math/lib/boost_1.78.0/boost/assert/source_location.hpp:80:18: warning: 'sprintf' is deprecated: This function is provided for compatibility reasons only.  Due to security concerns inherent in the design of sprintf(3), it is highly recommended that you use snprintf(3) instead. [-Wdeprecated-declarations]
            std::sprintf( buffer, ":%ld", static_cast<long>( column() ) );
                 ^
/Applications/Xcode.app/Contents/Developer/Platforms/MacOSX.platform/Developer/SDKs/MacOSX.sdk/usr/include/stdio.h:188:1: note: 'sprintf' has been explicitly marked deprecated here
__deprecated_msg("This function is provided for compatibility reasons only.  Due to security concerns inherent in the design of sprintf(3), it is highly recommended that you use snprintf(3) instead.")
^
/Applications/Xcode.app/Contents/Developer/Platforms/MacOSX.platform/Developer/SDKs/MacOSX.sdk/usr/include/sys/cdefs.h:215:48: note: expanded from macro '__deprecated_msg'
        #define __deprecated_msg(_msg) __attribute__((__deprecated__(_msg)))
                                                      ^
In file included from src/cmdstan/main.cpp:1:
In file included from src/cmdstan/command.hpp:4:
In file included from src/cmdstan/arguments/arg_data.hpp:4:
In file included from src/cmdstan/arguments/categorical_argument.hpp:4:
In file included from src/cmdstan/arguments/argument.hpp:4:
In file included from stan/src/stan/callbacks/writer.hpp:4:
In file included from stan/lib/stan_math/lib/boost_1.78.0/boost/lexical_cast.hpp:32:
In file included from stan/lib/stan_math/lib/boost_1.78.0/boost/lexical_cast/try_lexical_convert.hpp:44:
In file included from stan/lib/stan_math/lib/boost_1.78.0/boost/lexical_cast/detail/converter_lexical.hpp:54:
stan/lib/stan_math/lib/boost_1.78.0/boost/lexical_cast/detail/converter_lexical_streams.hpp:285:21: warning: 'sprintf' is deprecated: This function is provided for compatibility reasons only.  Due to security concerns inherent in the design of sprintf(3), it is highly recommended that you use snprintf(3) instead. [-Wdeprecated-declarations]
                    sprintf(begin,
                    ^
/Applications/Xcode.app/Contents/Developer/Platforms/MacOSX.platform/Developer/SDKs/MacOSX.sdk/usr/include/stdio.h:188:1: note: 'sprintf' has been explicitly marked deprecated here
__deprecated_msg("This function is provided for compatibility reasons only.  Due to security concerns inherent in the design of sprintf(3), it is highly recommended that you use snprintf(3) instead.")
^
/Applications/Xcode.app/Contents/Developer/Platforms/MacOSX.platform/Developer/SDKs/MacOSX.sdk/usr/include/sys/cdefs.h:215:48: note: expanded from macro '__deprecated_msg'
        #define __deprecated_msg(_msg) __attribute__((__deprecated__(_msg)))
                                                      ^
In file included from src/cmdstan/main.cpp:1:
In file included from src/cmdstan/command.hpp:4:
In file included from src/cmdstan/arguments/arg_data.hpp:4:
In file included from src/cmdstan/arguments/categorical_argument.hpp:4:
In file included from src/cmdstan/arg
uments/argument.hpp:4:
In file included from stan/src/stan/callbacks/writer.hpp:4:
In file included from stan/lib/stan_math/lib/boost_1.78.0/boost/lexical_cast.hpp:32:
In file included from stan/lib/stan_math/lib/boost_1.78.0/boost/lexical_cast/try_lexical_convert.hpp:44:
In file included from stan/lib/stan_math/lib/boost_1.78.0/boost/lexical_cast/detail/converter_lexical.hpp:54:
stan/lib/stan_math/lib/boost_1.78.0/boost/lexical_cast/detail/converter_lexical_streams.hpp:297:21: warning: 'sprintf' is deprecated: This function is provided for compatibility reasons only.  Due to security concerns inherent in the design of sprintf(3), it is highly recommended that you use snprintf(3) instead. [-Wdeprecated-declarations]
                    sprintf(begin,
                    ^
/Applications/Xcode.app/Contents/Developer/Platforms/MacOSX.platform/Developer/SDKs/MacOSX.sdk/usr/include/stdio.h:188:1: note: 'sprintf' has been explicitly marked deprecated here
__deprecated_msg("This function is provided for compatibility reasons only.  Due to security concerns inherent in the design of sprintf(3), it is highly recommended that you use snprintf(3) instead.")
^
/Applications/Xcode.app/Contents/Developer/Platforms/MacOSX.platform/Developer/SDKs/MacOSX.sdk/usr/include/sys/cdefs.h:215:48: note: 
expanded from macro '__deprecated_msg'
        #define __deprecated_msg(_msg) __attribute__((__deprecated__(_msg)))
                                                      ^
In file included from src/cmdstan/main.cpp:1:
In file included from src/cmdstan/command.hpp:4:
In file included from src/cmdstan/arguments/arg_data.hpp:4:
In file included from src/cmdstan/arguments/categorical_argument.hpp:4:
In file included from src/cmdstan/arguments/argument.hpp:4:
In file included from stan/src/stan/callbacks/writer.hpp:4:
In file included from stan/lib/stan_math/lib/boost_1.78.0/boost/lexical_cast.hpp:32:
In file included from stan/lib/stan_math/lib/boost_1.78.0/boost/lexical_cast/try_lexical_convert.hpp:44:
In file included from stan/lib/stan_math/lib/boost_1.78.0/boost/lexical_cast/detail/converter_lexical.hpp:54:
stan/lib/stan_math/lib/boost_1.78.0/boost/lexical_cast/detail/converter_lexical_streams.hpp:310:21: warning: 'sprintf' is deprecated: This function is provided for compatibility reasons only.  Due to security concerns inherent in the design of sprintf(3), it is highly recommended that you use snprintf(3) instead. [-Wdeprecated-declarations]
                    sprintf(begin,
                    ^
/Applications/Xcode.app/Contents/Developer/Platforms/MacOSX.platform/Developer/SDKs/MacOSX.sdk/usr/include/stdio.h:188:1: note: 'sprintf' has been explicitly marked deprecated here
__deprecated_msg("This function is provided for compatibility reasons only.  Due to security concerns inherent in the design of sprintf(3), it is highly recommended that you use snprintf(3) instead.")
^
/Applications/Xcode.app/Contents/Developer/Platforms/MacOSX.platform/Developer/SDKs/MacOSX.sdk/usr/include/sys/cdefs.h:215:48: note: expanded from macro '__deprecated_msg'
        #define __deprecated_msg(_msg) __attribute__((__deprecated__(_msg)))
                                                      ^
In file included from src/cmdstan/main.cpp:1:
In file included from src/cmdstan/command.hpp:29:
In file included from stan/src/stan/io/dump.hpp:7:
In file included from stan/lib/stan_math/stan/math/prim.hpp:15:
In file included from stan/lib/stan_math/stan/math/prim/functor.hpp:14:
In file included from stan/lib/stan_math/stan/math/prim/functor/integrate_ode_rk45.hpp:6:
In file included from stan/lib/stan_math/stan/math/prim/functor/ode_rk45.hpp:9:
In file included from stan/lib/stan_math/lib/boost_1.78.0/boost/numeric/odeint.hpp:37:
In file included from stan/lib/stan_math/lib/boost_1.78.0/boost/numeric/odeint/stepper/dense_output_runge_kutta.hpp:40:
stan/lib/stan_math/lib/boost_1.78.0/boost/numeric/odeint/integrate/max_step_checker.hpp:72:18: warning: 'sprintf' is deprecated: This function is provided for compatibility reasons only.  Due to security concerns inherent in the design of sprintf(3), it is highly recommended that you use snprintf(3) instead. [-Wdeprecated-declarations]
            std::sprintf(error_msg, "Max number of iterations exceeded (%d).", m_max_steps);
                 ^
/Applications/Xcode.app/Contents/Developer/Platforms/MacOSX.platform/Developer/SDKs/MacOSX.sdk/usr/include/stdio.h:188:1: note: 'sprintf' has been explicitly marked deprecated here
__deprecated_msg("This function is provided for compatibility reasons only.  Due to security concerns inherent in the design of sprintf(3), it is highly recommended that you use snprintf(3) instead.")
^
/Applications/Xcode.app/Contents/Developer/Platforms/MacOSX.platform/Developer/SDKs/MacOSX.sdk/usr/include/sys/cdefs.h:215:48: note: expanded from macro '__deprecated_msg'
        #define __deprecated_msg(_msg) __attribute__((__deprecated__(_msg)))
                                                      ^
In file included from src/cmdstan/main.cpp:1:
In file included from src/cmdstan/command.hpp:29:
In file included from stan/src/stan/io/dump.hpp:7:
In file included from stan/lib/stan_math/stan/math/prim.hpp:15:
In file included from stan/lib/stan_math/stan/math/prim/functor.hpp:14:
In file included from stan/lib/stan_math/stan/math/prim/functor/integrate_ode_rk45.hpp:6:
In file included from stan/lib/stan_math/stan/math/prim/functor/ode_rk45.hpp:9:
In file included from stan/lib/stan_math/lib/boost_1.78.0/boost/numeric/odeint.hpp:37:
In file included from stan/lib/stan_math/lib/boost_1.78.0/boost/numeric/odeint/stepper/dense_output_runge_kutta.hpp:40:
stan/lib/stan_math/lib/boost_1.78.0/boost/numeric/odeint/integrate/max_step_checker.hpp:104:18: warning: 'sprintf' is deprecated: This function is provided for compatibility reasons only.  Due to security concerns inherent in the design of sprintf(3), it is highly recommended that you use snprintf(3) instead. [-Wdeprecated-declarations]
            std::sprintf(error_msg, "Max number of iterations exceeded (%d). A new step size was not found.", m_max_steps);
                 ^
/Applications/Xcode.app/Contents/Developer/Platforms/MacOSX.platform/Developer/SDKs/MacOSX.sdk/usr/include/stdio.h:188:1: note: 'sprintf' has been explicitly marked deprecated here
__deprecated_msg("This function is provided for compatibility reasons only.  Due to security concerns inherent in the design of sprintf(3), it is highly recommended that you use snprintf(3) instead.")
^
/Applications/Xcode.app/Contents/Developer/Platforms/MacOSX.platform/Developer/SDKs/MacOSX.sdk/usr/include/sys/cdefs.h:215:48: note: expanded from macro '__deprecated_msg'
        #define __deprecated_msg(_msg) __attribute__((__deprecated__(_msg)))
                                                      ^
7 warnings generated.
In file included from stan/src/stan/model/model_header.hpp:4:
In file included from stan/lib/stan_math/stan/math.hpp:19:
In file included from stan/lib/stan_math/stan/math/rev.hpp:8:
In file included from stan/lib/stan_math/stan/math/rev/core.hpp:10:
In file included from stan/lib/stan_math/stan/math/rev/core/chainable_object.hpp:6:
In file included from stan/lib/stan_math/stan/math/rev/core/typedefs.hpp:7:
In file included from stan/lib/stan_math/stan/math/rev/core/Eigen_NumTraits.hpp:5:
In file included from stan/lib/stan_math/stan/math/prim/core.hpp:4:
In file included from stan/lib/stan_math/stan/math/prim/core/init_threadpool_tbb.hpp:6:
In file included from stan/lib/stan_math/lib/boost_1.78.0/boost/lexical_cast.hpp:31:
In file included from stan/lib/stan_math/lib/boost_1.78.0/boost/lexical_cast/bad_lexical_cast.hpp:28:
In file included from stan/lib/stan_math/lib/boost_1.78.0/boost/throw_exception.hpp:24:
stan/lib/stan_math/lib/boost_1.78.0/boost/assert/source_location.hpp:75:14: warning: 'sprintf' is deprecated: This function is provided for compatibility reasons only.  Due to security concerns inherent in the design of sprintf(3), it is highly recommended that you use snprintf(3) instead. [-Wdeprecated-declarations]
        std::sprintf( buffer, ":%ld", static_cast<long>( line() ) );
             ^
/Applications/Xcode.app/Contents/Developer/Platforms/MacOSX.platform/Developer/SDKs/MacOSX.sdk/usr/include/stdio.h:188:1: note: 'sprintf' has been explicitly marked deprecated here
__deprecated_msg("This function is provided for compatibility reasons only.  Due to security concerns inherent in the design of sprintf(3), it is highly recommended that you use snprintf(3) instead.")
^
/Applications/Xcode.app/Contents/Developer/Platforms/MacOSX.platform/Developer/SDKs/MacOSX.sdk/usr/include/sys/cdefs.h:215:48: note: expanded from macro '__deprecated_msg'
        #define __deprecated_msg(_msg) __attribute__((__deprecated__(_msg)))
                                       
               ^
In file included from stan/src/stan/model/model_header.hpp:4:
In file included from stan/lib/stan_math/stan/math.hpp:19:
In file included from stan/lib/stan_math/stan/math/rev.hpp:8:
In file included from stan/lib/stan_math/stan/math/rev/core.hpp:10:
In file included from stan/lib/stan_math/stan/math/rev/core/chainable_object.hpp:6:
In file included from stan/lib/stan_math/stan/math/rev/core/typedefs.hpp:7:
In file included from stan/lib/stan_math/stan/math/rev/core/Eigen_NumTraits.hpp:5:
In file included from stan/lib/stan_math/stan/math/prim/core.hpp:4:
In file included from stan/lib/stan_math/stan/math/prim/core/init_threadpool_tbb.hpp:6:
In file included from stan/lib/stan_math/lib/boost_1.78.0/boost/lexical_cast.hpp:31:
In file included from stan/lib/stan_math/lib/boost_1.78.0/boost/lexical_cast/bad_lexical_cast.hpp:28:
In file included from stan/lib/stan_math/lib/boost_1.78.0/boost/throw_exception.hpp:24:
stan/lib/stan_math/lib/boost_1.78.0/boost/assert/source_location.hpp:80:18: warning: 'sprintf' is deprecated: This function is provided for compatibility reasons only.  Due to security concerns inherent in the design of sprintf(3), it is highly recommended that you use snprintf(3) instead. [-Wdeprecated-declarations]
            std::sprintf( buffer, ":%ld", static_cast<long>( column() ) );
                 ^
/Applications/Xcode.app/Contents/Developer/Platforms/MacOSX.platform/Developer/SDKs/MacOSX.sdk/usr/include/stdio.h:188:1: note: 'sprintf' has been explicitly marked deprecated here
__deprecated_msg("This function is provided for compatibility reasons only.  Due to security concerns inherent in the design of sprintf(3), it is highly recommended that you use snprintf(3) instead.")
^
/Applications/Xcode.app/Contents/Developer/Platforms/MacOSX.platform/Developer/SDKs/MacOSX.sdk/usr/include/sys/cdefs.h:215:48: note: expanded from macro '__deprecated_msg'
        #define __deprecated_msg(_msg) __attribute__((__deprecated__(_msg)))
                             
                         ^
In file included from stan/src/stan/model/model_header.hpp:4:
In file included from stan/lib/stan_math/stan/math.hpp:19:
In file included from stan/lib/stan_math/stan/math/rev.hpp:8:
In file included from stan/lib/stan_math/stan/math/rev/core.hpp:10:
In file included from stan/lib/stan_math/stan/math/rev/core/chainable_object.hpp:6:
In file included from stan/lib/stan_math/stan/math/rev/core/typedefs.hpp:7:
In file included from stan/lib/stan_math/stan/math/rev/core/Eigen_NumTraits.hpp:5:
In file included from stan/lib/stan_math/stan/math/prim/core.hpp:4:
In file included from stan/lib/stan_math/stan/math/prim/core/init_threadpool_tbb.hpp:6:
In file included from stan/lib/stan_math/lib/boost_1.78.0/boost/lexical_cast.hpp:32:
In file included from stan/lib/stan_math/lib/boost_1.78.0/boost/lexical_cast/try_lexical_convert.hpp:44:
In file included from stan/lib/stan_math/lib/boost_1.78.0/boost/lexical_cast/detail/converter_lexical.hpp:54:
stan/lib/stan_math/lib/boost_1.78.0/boost/lexical_cast/detail/converter_lexical_streams.hpp:285:21: warning: 'sprintf' is deprecated: This function is provided for compatibility reasons only.  Due to security concerns inherent in the design of sprintf(3), it is highly recommended that you use snprintf(3) instead. [-Wdeprecated-declarations]
                    sprintf(begin,
                    ^
/Applications/Xcode.app/Contents/Developer/Platforms/MacOSX.platform/Developer/SDKs/MacOSX.sdk/usr/include/stdio.h:188:1: note: 'sprintf' has been explicitly marked deprecated here
__deprecated_msg("This function is provided for compatibility reasons only.  Due to security concerns inherent in the design of sprintf(3), it is highly recommended that you use snprintf(3) instead.")
^
/Applications/Xcode.app/Contents/Developer/Platforms/MacOSX.platform/Developer/SDKs/MacOSX.sdk/usr/include/sys/cdefs.h:215:48: note: expanded from macro '__deprecated_msg'
        #define __deprecated_msg(_msg) __attribute__((__deprecated__(_msg)))
                 
                                     ^
In file included from stan/src/stan/model/model_header.hpp:4:
In file included from stan/lib/stan_math/stan/math.hpp:19:
In file included from stan/lib/stan_math/stan/math/rev.hpp:8:
In file included from stan/lib/stan_math/stan/math/rev/core.hpp:10:
In file included from stan/lib/stan_math/stan/math/rev/core/chainable_object.hpp:6:
In file included from stan/lib/stan_math/stan/math/rev/core/typedefs.hpp:7:
In file included from stan/lib/stan_math/stan/math/rev/core/Eigen_NumTraits.hpp:5:
In file included from stan/lib/stan_math/stan/math/prim/core.hpp:4:
In file included from stan/lib/stan_math/stan/math/prim/core/init_threadpool_tbb.hpp:6:
In file included from stan/lib/stan_math/lib/boost_1.78.0/boost/lexical_cast.hpp:32:
In file included from stan/lib/stan_math/lib/boost_1.78.0/boost/lexical_cast/try_lexical_convert.hpp:44:
In file included from stan/lib/stan_math/lib/boost_1.78.0/boost/lexical_cast/detail/converter_lexical.hpp:54:
stan/lib/stan_math/lib/boost_1.78.0/boost/lexical_cast/detail/converter_lexical_streams.hpp:297:21: warning: 'sprintf' is deprecated: This function is provided for compatibility reasons only.  Due to security concerns inherent in the design of sprintf(3), it is highly recommended that you use snprintf(3) instead. [-Wdeprecated-declarations]
                    sprintf(begin,
                    ^
/Applications/Xcode.app/Contents/Developer/Platforms/MacOSX.platform/Developer/SDKs/MacOSX.sdk/usr/include/stdio.h:188:1: note: 'sprintf' has been explicitly marked deprecated here
__deprecated_msg("This function is provided for compatibility reasons only.  Due to security concerns inherent in the design of sprintf(3), it is highly recommended that you use snprintf(3) instead.")
^
/Applications/Xcode.app/Contents/Developer/Platforms/MacOSX.platform/Developer/SDKs/MacOSX.sdk/usr/include/sys/cdefs.h:215:48: note: expanded from macro '__deprecated_msg'
        #define __deprecated_msg(_msg) __attribute__((__deprecated__(_msg)))
                                                      ^
In file included from s
tan/src/stan/model/model_header.hpp:4:
In file included from stan/lib/stan_math/stan/math.hpp:19:
In file included from stan/lib/stan_math/stan/math/rev.hpp:8:
In file included from stan/lib/stan_math/stan/math/rev/core.hpp:10:
In file included from stan/lib/stan_math/stan/math/rev/core/chainable_object.hpp:6:
In file included from stan/lib/stan_math/stan/math/rev/core/typedefs.hpp:7:
In file included from stan/lib/stan_math/stan/math/rev/core/Eigen_NumTraits.hpp:5:
In file included from stan/lib/stan_math/stan/math/prim/core.hpp:4:
In file included from stan/lib/stan_math/stan/math/prim/core/init_threadpool_tbb.hpp:6:
In file included from stan/lib/stan_math/lib/boost_1.78.0/boost/lexical_cast.hpp:32:
In file included from stan/lib/stan_math/lib/boost_1.78.0/boost/lexical_cast/try_lexical_convert.hpp:44:
In file included from stan/lib/stan_math/lib/boost_1.78.0/boost/lexical_cast/detail/converter_lexical.hpp:54:
stan/lib/stan_math/lib/boost_1.78.0/boost/lexical_cast/detail/converter_lexical_streams.hpp:310:21: warning: 'sprintf' is deprecated: This function is provided for compatibility reasons only.  Due to security concerns inherent in the design of sprintf(3), it is highly recommended that you use snprintf(3) instead. [-Wdeprecated-declarations]
                    sprintf(begin,
                    ^
/Applications/Xcode.app/Contents/Developer/Platforms/MacOSX.platform/Developer/SDKs/MacOSX.sdk/usr/include/stdio.h:188:1: note: 'sprintf' has been explicitly marked deprecated here
__deprecated_msg("This function is provided for compatibility reasons only.  Due to security concerns inherent in the design of sprintf(3), it is highly recommended that you use snprintf(3) instead.")
^
/Applications/Xcode.app/Contents/Developer/Platforms/MacOSX.platform/Developer/SDKs/MacOSX.sdk/usr/include/sys/cdefs.h:215:48: note: expanded from macro '__deprecated_msg'
        #define __deprecated_msg(_msg) __attribute__((__deprecated__(_msg)))
                                        
              ^
In file included from stan/src/stan/model/model_header.hpp:4:
In file included from stan/lib/stan_math/stan/math.hpp:19:
In file included from stan/lib/stan_math/stan/math/rev.hpp:10:
In file included from stan/lib/stan_math/stan/math/rev/fun.hpp:198:
In file included from stan/lib/stan_math/stan/math/prim/functor.hpp:14:
In file included from stan/lib/stan_math/stan/math/prim/functor/integrate_ode_rk45.hpp:6:
In file included from stan/lib/stan_math/stan/math/prim/functor/ode_rk45.hpp:9:
In file included from stan/lib/stan_math/lib/boost_1.78.0/boost/numeric/odeint.hpp:37:
In file included from stan/lib/stan_math/lib/boost_1.78.0/boost/numeric/odeint/stepper/dense_output_runge_kutta.hpp:40:
stan/lib/stan_math/lib/boost_1.78.0/boost/numeric/odeint/integrate/max_step_checker.hpp:72:18: warning: 'sprintf' is deprecated: This function is provided for compatibility reasons only.  Due to security concerns inherent in the design of sprintf(3), it is highly recommended that you use snprintf(3) instead. [-Wdeprecated-declarations]
            std::sprintf(error_msg, "Max number of iterations exceeded (%d).", m_max_steps);
                 ^
/Applications/Xcode.app/Contents/Developer/Platforms/MacOSX.platform/Developer/SDKs/MacOSX.sdk/usr/include/stdio.h:188:1: note: 'sprintf' has been explicitly marked deprecated here
__deprecated_msg("This function is provided for compatibility reasons only.  Due to security concerns inherent in the design of sprintf(3), it is highly recommended that you use snprintf(3) instead.")
^
/Applications/Xcode.app/Contents/Developer/Platforms/MacOSX.platform/Developer/SDKs/MacOSX.sdk/usr/include/sys/cdefs.h:215:48: note: expanded from macro '__deprecated_msg'
        #define __deprecated_msg(_msg) __attribute__((__deprecated__(_msg)))
                                                      ^
In file included from stan/src/stan/model/model_header.hpp:4:
In file included from stan/lib/stan_math/stan/math.hpp:19:
In file included from stan/lib/stan_math/stan/math/rev.hpp:10:
In file included from stan/lib/stan_math/stan/math/rev/fun.hpp:198:
In file included from stan/lib/stan_math/stan/math/prim/functor.hpp:14:
In file included from stan/lib/stan_math/stan/math/prim/functor/integrate_ode_rk45.hpp:6:
In file included from stan/lib/stan_math/stan/math/prim/functor/ode_rk45.hpp:9:
In file included from stan/lib/stan_math/lib/boost_1.78.0/boost/numeric/odeint.hpp:37:
In file included from stan/lib/stan_math/lib/boost_1.78.0/boost/numeric/odeint/stepper/dense_output_runge_kutta.hpp:40:
stan/lib/stan_math/lib/boost_1.78.0/boost/numeric/odeint/integrate/max_step_checker.hpp:104:18: warning: 'sprintf' is deprecated: This function is provided for compatibility reasons only.  Due to security concerns inherent in the design of sprintf(3), it is highly recommended that you use snprintf(3) instead. [-Wdeprecated-declarations]
            std::sprintf(error_msg, "Max number of iterations exceeded (%d). A new step size was not found.", m_max_steps);
                 ^
/Applications/Xcode.app/Contents/Developer/Platforms/MacOSX.platform/Developer/SDKs/MacOSX.sdk/usr/include/stdio.h:188:1: note: 'sprintf' has been explicitly marked deprecated here
__deprecated_msg("This function is provided for compatibility reasons only.  Due to security concerns inherent in the design of sprintf(3), it is highly recommended that you use snprintf(3) instead.")
^
/Applications/Xcode.app/Contents/Developer/Platforms/MacOSX.platform/Developer/SDKs/MacOSX.sdk/usr/include/sys/cdefs.h:215:48: note: expanded from macro '__deprecated_msg'
        #define __deprecated_msg(_msg) __attribute__((__deprecated__(_msg)))
                                                      ^
7 warnings generated.
Running MCMC with 4 parallel chains, with 1 thread(s) per chain...

Chain 1 Iteration:   1 / 1000 [  0%]  (Warmup) 
Chain 1 Iteration: 100 / 1000 [ 10%]  (Warmup) 
Chain 1 Iteration: 200 / 1000 [ 20%]  (Warmup) 
Chain 1 Iteration: 300 / 1000 [ 30%]  (Warmup) 
Chain 1 Iteration: 400 / 1000 [ 40%]  (Warmup) 
Chain 2 Iteration:   1 / 1000 [  0%]  (Warmup) 
Chain 2 Iteration: 100 / 1000 [ 10%]  (Warmup) 
Chain 2 Iteration: 200 / 1000 [ 20%]  (Warmup) 
Chain 2 Iteration: 300 / 1000 [ 30%]  (Warmup) 
Chain 2 Iteration: 400 / 1000 [ 40%]  (Warmup) 
Chain 2 Iteration: 500 / 1000 [ 50%]  (Warmup) 
Chain 2 Iteration: 501 / 1000 [ 50%]  (Sampling) 
Chain 2 Iteration: 600 / 1000 [ 60%]  (Sampling) 
Chain 2 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:
Chain 2 Exception: normal_lpdf: Scale parameter is 0, but must be positive! (in '/var/folders/f9/mxzvfl4s7x7c5j52qhsgms7c0000gr/T/Rtmp0MtGZ7/model-12c0c5e24fe0c.stan', line 22, column 4 to column 34)
Chain 2 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,
Chain 2 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.
Chain 2 
Chain 3 Iteration:   1 / 1000 [  0%]  (Warmup) 
Chain 3 Iteration: 100 / 1000 [ 10%]  (Warmup) 
Chain 3 Iteration: 200 / 1000 [ 20%]  (Warmup) 
Chain 3 Iteration: 300 / 1000 [ 30%]  (Warmup) 
Chain 3 Iteration: 400 / 1000 [ 40%]  (Warmup) 
Chain 3 Iteration: 500 / 1000 [ 50%]  (Warmup) 
Chain 3 Iteration: 501 / 1000 [ 50%]  (Sampling) 
Chain 3 Iteration: 600 / 1000 [ 60%]  (Sampling) 
Chain 3 Iteration: 700 / 1000 [ 70%]  (Sampling) 
Chain 4 Iteration:   1 / 1000 [  0%]  (Warmup) 
Chain 4 Iteration: 100 / 1000 [ 10%]  (Warmup) 
Chain 4 Iteration: 200 / 1000 [ 20%]  (Warmup) 
Chain 4 Iteration: 300 / 1000 [ 30%]  (Warmup) 
Chain 4 Iteration: 400 / 1000 [ 40%]  (Warmup) 
Chain 4 Iteration: 500 / 1000 [ 50%]  (Warmup) 
Chain 4 Iteration: 501 / 1000 [ 50%]  (Sampling) 
Chain 4 Iteration: 600 / 1000 [ 60%]  (Sampling) 
Chain 4 Iteration: 700 / 1000 [ 70%]  (Sampling) 
Chain 1 Iteration: 500 / 1000 [ 50%]  (Warmup) 
Chain 1 Iteration: 501 / 1000 [ 50%]  (Sampling) 
Chain 1 Iteration: 600 / 1000 [ 60%]  (Sampling) 
Chain 1 Iteration: 700 / 1000 [ 70%]  (Sampling) 
Chain 1 Iteration: 800 / 1000 [ 80%]  (Sampling) 
Chain 1 Iteration: 900 / 1000 [ 90%]  (Sampling) 
Chain 1 Iteration: 1000 / 1000 [100%]  (Sampling) 
Chain 2 Iteration: 700 / 1000 [ 70%]  (Sampling) 
Chain 2 Iteration: 800 / 1000 [ 80%]  (Sampling) 
Chain 2 Iteration: 900 / 1000 [ 90%]  (Sampling) 
Chain 2 Iteration: 1000 / 1000 [100%]  (Sampling) 
Chain 3 Iteration: 800 / 1000 [ 80%]  (Sampling) 
Chain 3 Iteration: 900 / 1000 [ 90%]  (Sampling) 
Chain 3 Iteration: 1000 / 1000 [100%]  (Sampling) 
Chain 4 Iteration: 800 / 1000 [ 80%]  (Sampling) 
Chain 4 Iteration: 900 / 1000 [ 90%]  (Sampling) 
Chain 4 Iteration: 1000 / 1000 [100%]  (Sampling) 
Chain 1 finished in 0.3 seconds.
Chain 2 finished in 0.4 seconds.
Chain 3 finished in 0.3 seconds.
Chain 4 finished in 0.3 seconds.

All 4 chains finished successfully.
Mean chain execution time: 0.3 seconds.
Total execution time: 0.5 seconds.

Let’s see the Stan code!

stancode(fit)
data{
    vector[116] area;
    vector[116] weight;
    vector[116] groupsize;
    vector[116] avgfood;
}
parameters{
    real beta0;
    real beta_group;
    real beta_food;
    real<lower=0> sigma;
}
model{
    vector[116] mu;
    sigma ~ exponential( 1 );
    beta_food ~ normal( 0 , 0.5 );
    beta_group ~ normal( 0 , 0.5 );
    beta0 ~ normal( 0 , 0.2 );
    for ( i in 1:116 ) {
        mu[i] = beta0 + beta_food * avgfood[i] + beta_group * groupsize[i];
    }
    weight ~ normal( mu , sigma );
}
generated quantities{
    vector[116] log_lik;
    vector[116] mu;
    for ( i in 1:116 ) {
        mu[i] = beta0 + beta_food * avgfood[i] + beta_group * groupsize[i];
    }
    for ( i in 1:116 ) log_lik[i] = normal_lpdf( weight[i] | mu[i] , sigma );
}

Data Block

Specify the data and indices. Note that “declaration of arrays by placing brackets after a variable name is deprecated.”

data {
  int<lower = 1> N; // Number of observations.
  int<lower = 1> I; // Number of covariates.
  vector[N] y;      // Vector of observations.
  matrix[N, I] X;   // Matrix of covariates.
}

Parameters Block

Specify the parameters (i.e., everything that is unknown).

parameters {
  vector[I] beta;        // Vector of regression parameters.
  real<lower = 0> sigma; // Variance of the likelihood.
}

Model Block

Specify the joint model: the priors and the likelihood (i.e., observational model).

model {
  // Priors.
  beta ~ normal(0, 5);
  sigma ~ normal(0, 5);

  // Likelihood.
  for (n in 1:N) {
    y[n] ~ normal(X[n,] * beta, sigma);
  }
}

(Optional) Generated Quantities Block

Specify anything you want to have efficiently computed along with the posterior draws. The most common generated quantity is the log-likelihood, which we’ll need to compute WAIC and PSIS.

generated quantities {
  // Compute the log-likelihood.
  vector[N] log_lik;
  for (n in 1:N) {
    log_lik[n] = normal_lpdf(y[n] | X[n,] * beta, sigma);
  }
}

We can use the data and generated quantities blocks to simulate data (kind of like extract.prior()).

data {
  int<lower = 1> N;               // Number of observations.
  int<lower = 1> I;               // Number of covariates.
  matrix[N, I] X;                 // Matrix of covariates.
  real<lower = 0> sigma;          // Variance of the likelihood.
}

generated quantities {
  vector[N] y;                    // Vector of observations.
  vector[I] beta;                 // Vector of regression parameters.

  // Draw parameter values and generate data.
  beta = normal_rng(0, sigma);
  for (n in 1:N) {
    y[n] = normal_rng(X[n,] * beta, sigma);
  }
}

RStudio and Stan

There is built-in support for Stan in RStudio.

  • Specify a stan code cell (requires #| output.var: output-name).
  • Create a .stan file (always leave an empty line at the end).

For simple models, you might get away with stan code cells. If you do, you’ll need to tell the knitr engine to use CmdStanR instead of RStan when working with stan code cells:

library(cmdstanr)
check_cmdstan_toolchain(fix = TRUE, quiet = TRUE)
register_knitr_engine()

However…

Stan and CmdStanR

As your project’s models get more complicated, we typically don’t want to run them as part of rendering a Quarto document. Instead, save them as their own .stan files, run them from a separate R script, save the model output, and print the model and import the output as part of our Quarto document as needed. You may also need to keep your computer awake long enough to run the model.

See the CmdStanR Vignette for more details. Imagine we’re doing this in a separate R script from our Quarto document so it doesn’t run every time we render.

First we need to compile the model (if it isn’t already).

# Load packages.
library(rethinking)
library(cmdstanr)
library(posterior)
library(bayesplot)
library(tidyverse)

# Compile the model.
stan_reg <- cmdstan_model(here::here("Code", "regression.stan"))

In a departure from typical R syntax, we will use $ notation to call functions. It works like this: model_object$function(). This strange syntax allows for consistency across CmdStan’s variants, including CmdStanR and CmdStanPy.

# Print the models like with stancode().
stan_reg$print()
// Index values, observations, and covariates.
data {
  int<lower = 1> N; // Number of observations.
  int<lower = 1> I; // Number of covariates.
  vector[N] y;      // Vector of observations.
  matrix[N, I] X;   // Matrix of covariates.
}

// Parameters.
parameters {
  vector[I] beta;        // Vector of regression parameters.
  real<lower = 0> sigma; // Variance of the likelihood.
}

// Regression model.
model {
  // Priors.
  beta ~ normal(0, 5);
  sigma ~ normal(0, 5);

  // Likelihood.
  for (n in 1:N) {
    y[n] ~ normal(X[n,] * beta, sigma);
  }
}

// Compute the log-likelihood.
generated quantities {
  // Compute the log-likelihood.
  vector[N] log_lik;
  for (n in 1:N) {
    log_lik[n] = normal_lpdf(y[n] | X[n,] * beta, sigma);
  }
}
# Structure the data to match the data block.
data_list <- list(
  N = nrow(foxes),
  I = 3,
  y = foxes$weight,
  X = as.matrix(bind_cols(1, foxes$avgfood, foxes$groupsize))
)
New names:
• `` -> `...1`
• `` -> `...2`
• `` -> `...3`
# Fit the model by sampling from the posterior.
fit <- stan_reg$sample(
  data = data_list, 
  seed = 123, 
  chains = 4, 
  parallel_chains = 4
)
Running MCMC with 4 parallel chains...

Chain 1 Iteration:    1 / 2000 [  0%]  (Warmup) 
Chain 1 Iteration:  100 / 2000 [  5%]  (Warmup) 
Chain 1 Iteration:  200 / 2000 [ 10%]  (Warmup) 
Chain 1 Iteration:  300 / 2000 [ 15%]  (Warmup) 
Chain 1 Iteration:  400 / 2000 [ 20%]  (Warmup) 
Chain 1 Iteration:  500 / 2000 [ 25%]  (Warmup) 
Chain 1 Iteration:  600 / 2000 [ 30%]  (Warmup) 
Chain 1 Iteration:  700 / 2000 [ 35%]  (Warmup) 
Chain 1 Iteration:  800 / 2000 [ 40%]  (Warmup) 
Chain 1 Iteration:  900 / 2000 [ 45%]  (Warmup) 
Chain 1 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:
Chain 1 Exception: normal_lpdf: Scale parameter is 0, but must be positive! (in '/var/folders/f9/mxzvfl4s7x7c5j52qhsgms7c0000gr/T/Rtmp5JnlM0/model-12890102c8b07.stan', line 23, column 4 to column 39)
Chain 1 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,
Chain 1 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.
Chain 1 
Chain 2 Iteration:    1 / 2000 [  0%]  (Warmup) 
Chain 2 Iteration:  100 / 2000 [  5%]  (Warmup) 
Chain 2 Iteration:  200 / 2000 [ 10%]  (Warmup) 
Chain 2 Iteration:  300 / 2000 [ 15%]  (Warmup) 
Chain 2 Iteration:  400 / 2000 [ 20%]  (Warmup) 
Chain 2 Iteration:  500 / 2000 [ 25%]  (Warmup) 
Chain 2 Iteration:  600 / 2000 [ 30%]  (Warmup) 
Chain 2 Iteration:  700 / 2000 [ 35%]  (Warmup) 
Chain 2 Iteration:  800 / 2000 [ 40%]  (Warmup) 
Chain 2 Iteration:  900 / 2000 [ 45%]  (Warmup) 
Chain 2 Iteration: 1000 / 2000 [ 50%]  (Warmup) 
Chain 2 Iteration: 1001 / 2000 [ 50%]  (Sampling) 
Chain 2 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:
Chain 2 Exception: normal_lpdf: Scale parameter is 0, but must be positive! (in '/var/folders/f9/mxzvfl4s7x7c5j52qhsgms7c0000gr/T/Rtmp5JnlM0/model-12890102c8b07.stan', line 23, column 4 to column 39)
Chain 2 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,
Chain 2 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.
Chain 2 
Chain 3 Iteration:    1 / 2000 [  0%]  (Warmup) 
Chain 3 Iteration:  100 / 2000 [  5%]  (Warmup) 
Chain 3 Iteration:  200 / 2000 [ 10%]  (Warmup) 
Chain 3 Iteration:  300 / 2000 [ 15%]  (Warmup) 
Chain 3 Iteration:  400 / 2000 [ 20%]  (Warmup) 
Chain 3 Iteration:  500 / 2000 [ 25%]  (Warmup) 
Chain 3 Iteration:  600 / 2000 [ 30%]  (Warmup) 
Chain 3 Iteration:  700 / 2000 [ 35%]  (Warmup) 
Chain 3 Iteration:  800 / 2000 [ 40%]  (Warmup) 
Chain 3 Iteration:  900 / 2000 [ 45%]  (Warmup) 
Chain 4 Iteration:    1 / 2000 [  0%]  (Warmup) 
Chain 4 Iteration:  100 / 2000 [  5%]  (Warmup) 
Chain 4 Iteration:  200 / 2000 [ 10%]  (Warmup) 
Chain 4 Iteration:  300 / 2000 [ 15%]  (Warmup) 
Chain 4 Iteration:  400 / 2000 [ 20%]  (Warmup) 
Chain 4 Iteration:  500 / 2000 [ 25%]  (Warmup) 
Chain 4 Iteration:  600 / 2000 [ 30%]  (Warmup) 
Chain 4 Iteration:  700 / 2000 [ 35%]  (Warmup) 
Chain 4 Iteration:  800 / 2000 [ 40%]  (Warmup) 
Chain 4 Iteration:  900 / 2000 [ 45%]  (Warmup) 
Chain 4 Iteration: 1000 / 2000 [ 50%]  (Warmup) 
Chain 4 Iteration: 1001 / 2000 [ 50%]  (Sampling) 
Chain 1 Iteration: 1000 / 2000 [ 50%]  (Warmup) 
Chain 1 Iteration: 1001 / 2000 [ 50%]  (Sampling) 
Chain 1 Iteration: 1100 / 2000 [ 55%]  (Sampling) 
Chain 1 Iteration: 1200 / 2000 [ 60%]  (Sampling) 
Chain 1 Iteration: 1300 / 2000 [ 65%]  (Sampling) 
Chain 1 Iteration: 1400 / 2000 [ 70%]  (Sampling) 
Chain 1 Iteration: 1500 / 2000 [ 75%]  (Sampling) 
Chain 1 Iteration: 1600 / 2000 [ 80%]  (Sampling) 
Chain 2 Iteration: 1100 / 2000 [ 55%]  (Sampling) 
Chain 2 Iteration: 1200 / 2000 [ 60%]  (Sampling) 
Chain 2 Iteration: 1300 / 2000 [ 65%]  (Sampling) 
Chain 2 Iteration: 1400 / 2000 [ 70%]  (Sampling) 
Chain 2 Iteration: 1500 / 2000 [ 75%]  (Sampling) 
Chain 2 Iteration: 1600 / 2000 [ 80%]  (Sampling) 
Chain 3 Iteration: 1000 / 2000 [ 50%]  (Warmup) 
Chain 3 Iteration: 1001 / 2000 [ 50%]  (Sampling) 
Chain 3 Iteration: 1100 / 2000 [ 55%]  (Sampling) 
Chain 3 Iteration: 1200 / 2000 [ 60%]  (Sampling) 
Chain 3 Iteration: 1300 / 2000 [ 65%]  (Sampling) 
Chain 3 Iteration: 1400 / 2000 [ 70%]  (Sampling) 
Chain 3 Iteration: 1500 / 2000 [ 75%]  (Sampling) 
Chain 4 Iteration: 1100 / 2000 [ 55%]  (Sampling) 
Chain 4 Iteration: 1200 / 2000 [ 60%]  (Sampling) 
Chain 4 Iteration: 1300 / 2000 [ 65%]  (Sampling) 
Chain 4 Iteration: 1400 / 2000 [ 70%]  (Sampling) 
Chain 4 Iteration: 1500 / 2000 [ 75%]  (Sampling) 
Chain 1 Iteration: 1700 / 2000 [ 85%]  (Sampling) 
Chain 1 Iteration: 1800 / 2000 [ 90%]  (Sampling) 
Chain 1 Iteration: 1900 / 2000 [ 95%]  (Sampling) 
Chain 1 Iteration: 2000 / 2000 [100%]  (Sampling) 
Chain 2 Iteration: 1700 / 2000 [ 85%]  (Sampling) 
Chain 2 Iteration: 1800 / 2000 [ 90%]  (Sampling) 
Chain 2 Iteration: 1900 / 2000 [ 95%]  (Sampling) 
Chain 2 Iteration: 2000 / 2000 [100%]  (Sampling) 
Chain 3 Iteration: 1600 / 2000 [ 80%]  (Sampling) 
Chain 3 Iteration: 1700 / 2000 [ 85%]  (Sampling) 
Chain 3 Iteration: 1800 / 2000 [ 90%]  (Sampling) 
Chain 3 Iteration: 1900 / 2000 [ 95%]  (Sampling) 
Chain 3 Iteration: 2000 / 2000 [100%]  (Sampling) 
Chain 4 Iteration: 1600 / 2000 [ 80%]  (Sampling) 
Chain 4 Iteration: 1700 / 2000 [ 85%]  (Sampling) 
Chain 4 Iteration: 1800 / 2000 [ 90%]  (Sampling) 
Chain 4 Iteration: 1900 / 2000 [ 95%]  (Sampling) 
Chain 4 Iteration: 2000 / 2000 [100%]  (Sampling) 
Chain 1 finished in 0.3 seconds.
Chain 2 finished in 0.3 seconds.
Chain 3 finished in 0.4 seconds.
Chain 4 finished in 0.3 seconds.

All 4 chains finished successfully.
Mean chain execution time: 0.3 seconds.
Total execution time: 0.4 seconds.
# Look at model diagnostics.
fit$cmdstan_diagnose()
Processing csv files: /var/folders/f9/mxzvfl4s7x7c5j52qhsgms7c0000gr/T/Rtmp0MtGZ7/regression-202303151112-1-6484b3.csv, /var/folders/f9/mxzvfl4s7x7c5j52qhsgms7c0000gr/T/Rtmp0MtGZ7/regression-202303151112-2-6484b3.csv, /var/folders/f9/mxzvfl4s7x7c5j52qhsgms7c0000gr/T/Rtmp0MtGZ7/regression-202303151112-3-6484b3.csv, /var/folders/f9/mxzvfl4s7x7c5j52qhsgms7c0000gr/T/Rtmp0MtGZ7/regression-202303151112-4-6484b3.csv

Checking sampler transitions treedepth.
Treedepth satisfactory for all transitions.

Checking sampler transitions for divergences.
No divergent transitions found.

Checking E-BFMI - sampler transitions HMC potential energy.
E-BFMI satisfactory.

Effective sample size satisfactory.

Split R-hat values satisfactory all parameters.

Processing complete, no problems detected.
# Summarize model fit, including referencing specific parameters, like precis().
fit$cmdstan_summary()
Inference for Stan model: regression_model
4 chains: each with iter=(1000,1000,1000,1000); warmup=(0,0,0,0); thin=(1,1,1,1); 4000 iterations saved.

Warmup took (0.13, 0.13, 0.14, 0.14) seconds, 0.53 seconds total
Sampling took (0.21, 0.22, 0.22, 0.21) seconds, 0.87 seconds total

                    Mean     MCSE  StdDev     5%       50%    95%    N_Eff  N_Eff/s    R_hat

lp__            -5.3e+01  3.4e-02     1.4    -56  -5.3e+01    -52     1665     1922      1.0
accept_stat__       0.92  1.7e-03    0.11   0.68      0.96    1.0  4.3e+03  5.0e+03  1.0e+00
stepsize__          0.33  1.5e-02   0.021   0.31      0.33   0.37  2.0e+00  2.3e+00  1.8e+13
treedepth__          3.2  2.1e-02    0.71    2.0       3.0    4.0  1.2e+03  1.3e+03  1.0e+00
n_leapfrog__          11  3.2e-01     4.4    3.0        15     15  1.8e+02  2.1e+02  1.0e+00
divergent__         0.00      nan    0.00   0.00      0.00   0.00      nan      nan      nan
energy__              55  5.1e-02     2.0     53        55     59  1.5e+03  1.8e+03  1.0e+00

beta[1]         -3.9e-04  1.6e-03   0.088  -0.15  -1.8e-03   0.15     3152     3640      1.0
beta[2]          6.4e-01  4.2e-03    0.20   0.31   6.4e-01   0.95     2371     2738      1.0
beta[3]         -7.4e-01  4.1e-03    0.20   -1.1  -7.4e-01  -0.40     2432     2809      1.0
sigma            9.6e-01  1.2e-03   0.064   0.86   9.6e-01    1.1     2824     3261     1.00
log_lik[1]      -1.0e+00  2.0e-03    0.12   -1.3  -1.0e+00  -0.87     3658     4224     1.00
log_lik[2]      -1.9e+00  4.6e-03    0.29   -2.3  -1.8e+00   -1.4     3943     4553     1.00
log_lik[3]      -9.3e-01  1.5e-03   0.081   -1.1  -9.3e-01  -0.81     2960     3418      1.0
log_lik[4]      -1.3e+00  2.9e-03    0.17   -1.6  -1.3e+00   -1.1     3416     3944      1.0
log_lik[5]      -1.3e+00  2.4e-03    0.15   -1.5  -1.3e+00   -1.1     3851     4447     1.00
log_lik[6]      -1.9e+00  3.9e-03    0.25   -2.3  -1.9e+00   -1.5     4002     4621     1.00
log_lik[7]      -9.1e-01  1.4e-03   0.074   -1.0  -9.0e-01  -0.79     2713     3133      1.0
log_lik[8]      -1.0e+00  1.8e-03    0.11   -1.2  -1.0e+00  -0.88     3539     4086      1.0
log_lik[9]      -1.2e+00  3.2e-03    0.16   -1.5  -1.2e+00  -0.98     2420     2794      1.0
log_lik[10]     -9.5e-01  1.8e-03   0.089   -1.1  -9.4e-01  -0.81     2482     2866      1.0
log_lik[11]     -1.0e+00  2.3e-03    0.12   -1.2  -1.0e+00  -0.87     2453     2833      1.0
log_lik[12]     -1.1e+00  1.5e-03   0.087   -1.2  -1.1e+00  -0.93     3286     3794      1.0
log_lik[13]     -9.1e-01  1.3e-03   0.071   -1.0  -9.1e-01  -0.80     2771     3199      1.0
log_lik[14]     -1.8e+00  2.9e-03    0.17   -2.1  -1.7e+00   -1.5     3512     4055     1.00
log_lik[15]     -1.9e+00  6.0e-03    0.29   -2.4  -1.9e+00   -1.5     2403     2775      1.0
log_lik[16]     -1.2e+00  3.0e-03    0.15   -1.4  -1.1e+00  -0.95     2427     2803      1.0
log_lik[17]     -2.5e+00  7.2e-03    0.38   -3.2  -2.5e+00   -1.9     2786     3217      1.0
log_lik[18]     -9.4e-01  1.8e-03   0.093   -1.1  -9.3e-01  -0.80     2733     3156      1.0
log_lik[19]     -9.1e-01  1.5e-03   0.077   -1.0  -9.1e-01  -0.79     2601     3003     1.00
log_lik[20]     -2.0e+00  6.1e-03    0.34   -2.6  -2.0e+00   -1.5     3022     3490      1.0
log_lik[21]     -1.2e+00  2.4e-03    0.13   -1.5  -1.2e+00   -1.0     2731     3153      1.0
log_lik[22]     -9.4e-01  1.6e-03   0.080   -1.1  -9.4e-01  -0.82     2537     2929      1.0
log_lik[23]     -1.5e+00  3.2e-03    0.16   -1.8  -1.4e+00   -1.2     2687     3102      1.0
log_lik[24]     -2.1e+00  4.5e-03    0.24   -2.6  -2.1e+00   -1.8     2854     3295      1.0
log_lik[25]     -1.6e+00  3.0e-03    0.16   -1.8  -1.5e+00   -1.3     2974     3435      1.0
log_lik[26]     -4.3e+00  9.3e-03    0.52   -5.3  -4.3e+00   -3.5     3170     3660      1.0
log_lik[27]     -8.9e-01  1.3e-03   0.069   -1.0  -8.9e-01  -0.78     2713     3133     1.00
log_lik[28]     -1.2e+00  2.2e-03    0.13   -1.5  -1.2e+00   -1.0     3503     4045     1.00
log_lik[29]     -1.9e+00  3.8e-03    0.22   -2.3  -1.9e+00   -1.5     3397     3923     1.00
log_lik[30]     -1.7e+00  4.3e-03    0.25   -2.1  -1.6e+00   -1.3     3488     4028      1.0
log_lik[31]     -9.0e-01  1.4e-03   0.073   -1.0  -9.0e-01  -0.79     2590     2991     1.00
log_lik[32]     -1.1e+00  2.3e-03    0.13   -1.3  -1.1e+00  -0.91     3357     3877      1.0
log_lik[33]     -1.4e+00  3.4e-03    0.20   -1.7  -1.4e+00   -1.1     3420     3949      1.0
log_lik[34]     -1.2e+00  2.8e-03    0.16   -1.5  -1.2e+00  -0.98     3420     3949      1.0
log_lik[35]     -9.4e-01  1.6e-03   0.087   -1.1  -9.3e-01  -0.81     2944     3399      1.0
log_lik[36]     -4.3e+00  1.1e-02    0.65   -5.5  -4.3e+00   -3.3     3321     3835      1.0
log_lik[37]     -8.9e-01  1.3e-03   0.069   -1.0  -8.9e-01  -0.78     2685     3100     1.00
log_lik[38]     -1.1e+00  1.9e-03    0.10   -1.2  -1.1e+00  -0.92     2854     3296      1.0
log_lik[39]     -8.9e-01  1.3e-03   0.068   -1.0  -8.9e-01  -0.78     2679     3094     1.00
log_lik[40]     -2.1e+00  4.8e-03    0.26   -2.6  -2.1e+00   -1.8     2838     3277     1.00
log_lik[41]     -8.9e-01  1.3e-03   0.068  -1.00  -8.9e-01  -0.78     2633     3041     1.00
log_lik[42]     -9.4e-01  1.4e-03   0.072   -1.1  -9.4e-01  -0.83     2587     2987      1.0
log_lik[43]     -8.9e-01  1.3e-03   0.067  -1.00  -8.9e-01  -0.78     2673     3087     1.00
log_lik[44]     -1.5e+00  2.4e-03    0.12   -1.7  -1.5e+00   -1.3     2690     3106      1.0
log_lik[45]     -1.4e+00  1.8e-03    0.10   -1.6  -1.4e+00   -1.3     3060     3534      1.0
log_lik[46]     -9.4e-01  1.2e-03   0.067   -1.1  -9.4e-01  -0.83     2908     3358     1.00
log_lik[47]     -1.1e+00  1.3e-03   0.072   -1.2  -1.1e+00  -0.98     3100     3580      1.0
log_lik[48]     -3.6e+00  6.5e-03    0.36   -4.2  -3.6e+00   -3.0     3092     3571     1.00
log_lik[49]     -9.3e-01  1.3e-03   0.068   -1.0  -9.3e-01  -0.82     2712     3131     1.00
log_lik[50]     -9.0e-01  1.3e-03   0.068   -1.0  -9.0e-01  -0.79     2697     3114     1.00
log_lik[51]     -9.1e-01  1.2e-03   0.066   -1.0  -9.1e-01  -0.80     2829     3267     1.00
log_lik[52]     -1.7e+00  2.1e-03    0.12   -1.9  -1.7e+00   -1.5     3317     3830      1.0
log_lik[53]     -9.1e-01  1.3e-03   0.068   -1.0  -9.1e-01  -0.80     2823     3260     1.00
log_lik[54]     -9.9e-01  1.4e-03   0.073   -1.1  -9.9e-01  -0.87     2880     3326     1.00
log_lik[55]     -9.6e-01  1.3e-03   0.070   -1.1  -9.6e-01  -0.85     2926     3379     1.00
log_lik[56]     -2.7e+00  5.1e-03    0.27   -3.2  -2.7e+00   -2.3     2916     3368      1.0
log_lik[57]     -1.2e+00  3.0e-03    0.20   -1.6  -1.2e+00  -0.96     4234     4889      1.0
log_lik[58]     -9.3e-01  1.7e-03   0.089   -1.1  -9.2e-01  -0.80     2809     3243     1.00
log_lik[59]     -9.6e-01  1.9e-03    0.11   -1.2  -9.4e-01  -0.82     3250     3753     1.00
log_lik[60]     -9.3e-01  1.7e-03   0.090   -1.1  -9.2e-01  -0.80     2836     3274     1.00
log_lik[61]     -9.3e-01  1.7e-03   0.091   -1.1  -9.2e-01  -0.80     2845     3286      1.0
log_lik[62]     -1.0e+00  2.1e-03    0.13   -1.3  -1.0e+00  -0.85     3724     4300      1.0
log_lik[63]     -1.0e+00  2.1e-03    0.13   -1.2  -9.9e-01  -0.84     3643     4206     1.00
log_lik[64]     -1.0e+00  2.2e-03    0.13   -1.3  -1.0e+00  -0.85     3750     4330      1.0
log_lik[65]     -1.1e+00  1.3e-03   0.075   -1.2  -1.1e+00  -0.99     3178     3670     1.00
log_lik[66]     -1.6e+00  2.1e-03    0.12   -1.8  -1.6e+00   -1.4     3155     3643     1.00
log_lik[67]     -8.9e-01  1.3e-03   0.067  -1.00  -8.9e-01  -0.78     2757     3183     1.00
log_lik[68]     -1.4e+00  1.8e-03    0.10   -1.6  -1.4e+00   -1.3     3333     3849     1.00
log_lik[69]     -1.4e+00  2.3e-03    0.12   -1.6  -1.4e+00   -1.2     2748     3173      1.0
log_lik[70]     -1.4e+00  2.2e-03    0.12   -1.6  -1.4e+00   -1.2     2736     3159      1.0
log_lik[71]     -1.0e+00  1.6e-03   0.079   -1.1  -1.0e+00  -0.88     2578     2977      1.0
log_lik[72]     -1.5e+00  2.5e-03    0.13   -1.7  -1.5e+00   -1.3     2783     3214      1.0
log_lik[73]     -1.3e+00  2.2e-03    0.11   -1.5  -1.3e+00   -1.2     2728     3150      1.0
log_lik[74]     -2.7e+00  5.2e-03    0.28   -3.2  -2.7e+00   -2.3     2982     3443      1.0
log_lik[75]     -1.0e+00  1.5e-03   0.079   -1.2  -1.0e+00  -0.91     2768     3196      1.0
log_lik[76]     -1.7e+00  3.2e-03    0.17   -2.0  -1.7e+00   -1.5     2588     2989      1.0
log_lik[77]     -1.2e+00  2.1e-03    0.11   -1.4  -1.2e+00   -1.0     2572     2970      1.0
log_lik[78]     -9.2e-01  1.4e-03   0.073   -1.0  -9.2e-01  -0.80     2517     2906     1.00
log_lik[79]     -9.6e-01  1.5e-03   0.076   -1.1  -9.6e-01  -0.84     2712     3131      1.0
log_lik[80]     -3.6e+00  8.6e-03    0.44   -4.4  -3.6e+00   -2.9     2608     3011      1.0
log_lik[81]     -9.2e-01  1.2e-03   0.067   -1.0  -9.1e-01  -0.81     2859     3302     1.00
log_lik[82]     -1.3e+00  1.5e-03   0.087   -1.4  -1.3e+00   -1.1     3266     3772     1.00
log_lik[83]     -8.9e-01  1.3e-03   0.067   -1.0  -8.9e-01  -0.78     2730     3153     1.00
log_lik[84]     -1.0e+00  1.2e-03   0.068   -1.1  -1.0e+00  -0.89     3049     3521     1.00
log_lik[85]     -9.7e-01  1.3e-03   0.069   -1.1  -9.7e-01  -0.86     2779     3209      1.0
log_lik[86]     -9.7e-01  1.3e-03   0.069   -1.1  -9.6e-01  -0.86     2770     3198      1.0
log_lik[87]     -1.2e+00  1.4e-03   0.078   -1.3  -1.2e+00   -1.1     3240     3741      1.0
log_lik[88]     -2.0e+00  3.0e-03    0.16   -2.3  -2.0e+00   -1.7     2961     3419      1.0
log_lik[89]     -9.9e-01  2.4e-03    0.12   -1.2  -9.8e-01  -0.83     2446     2824     1.00
log_lik[90]     -9.1e-01  1.7e-03   0.081   -1.0  -9.1e-01  -0.79     2352     2716      1.0
log_lik[91]     -9.2e-01  1.7e-03   0.085   -1.1  -9.1e-01  -0.79     2397     2768      1.0
log_lik[92]     -1.6e+00  5.7e-03    0.29   -2.2  -1.6e+00   -1.2     2564     2960      1.0
log_lik[93]     -2.6e+00  8.8e-03    0.45   -3.4  -2.5e+00   -1.9     2652     3062      1.0
log_lik[94]     -1.2e+00  1.7e-03   0.092   -1.4  -1.2e+00   -1.1     2873     3317      1.0
log_lik[95]     -9.0e-01  1.3e-03   0.067   -1.0  -9.0e-01  -0.79     2733     3156     1.00
log_lik[96]     -9.4e-01  1.3e-03   0.070   -1.1  -9.4e-01  -0.83     2772     3201     1.00
log_lik[97]     -9.7e-01  1.3e-03   0.072   -1.1  -9.7e-01  -0.85     2873     3317     1.00
log_lik[98]     -1.2e+00  1.6e-03   0.088   -1.3  -1.2e+00   -1.0     3078     3554      1.0
log_lik[99]     -9.0e-01  1.3e-03   0.068   -1.0  -9.0e-01  -0.79     2753     3179     1.00
log_lik[100]    -1.3e+00  1.7e-03   0.097   -1.5  -1.3e+00   -1.1     3198     3693     1.00
log_lik[101]    -9.4e-01  1.3e-03   0.069   -1.1  -9.4e-01  -0.83     2830     3268     1.00
log_lik[102]    -1.6e+00  2.3e-03    0.13   -1.8  -1.6e+00   -1.4     3175     3667     1.00
log_lik[103]    -2.1e+00  3.4e-03    0.19   -2.4  -2.1e+00   -1.8     3044     3515     1.00
log_lik[104]    -1.1e+00  1.8e-03    0.10   -1.3  -1.1e+00  -0.98     3296     3806      1.0
log_lik[105]    -9.0e-01  1.3e-03   0.070   -1.0  -9.0e-01  -0.79     2717     3137     1.00
log_lik[106]    -4.2e+00  9.1e-03    0.52   -5.1  -4.1e+00   -3.4     3182     3675      1.0
log_lik[107]    -1.9e+00  3.6e-03    0.22   -2.3  -1.9e+00   -1.6     3531     4077      1.0
log_lik[108]    -1.7e+00  3.2e-03    0.19   -2.1  -1.7e+00   -1.5     3572     4125      1.0
log_lik[109]    -2.3e+00  4.5e-03    0.26   -2.7  -2.2e+00   -1.9     3254     3758      1.0
log_lik[110]    -1.0e+00  1.3e-03   0.072   -1.1  -1.0e+00  -0.91     3032     3501     1.00
log_lik[111]    -1.0e+00  1.4e-03   0.074   -1.2  -1.0e+00  -0.91     2938     3392     1.00
log_lik[112]    -9.5e-01  1.3e-03   0.068   -1.1  -9.4e-01  -0.84     2917     3368     1.00
log_lik[113]    -9.7e-01  1.3e-03   0.071   -1.1  -9.7e-01  -0.86     2850     3291     1.00
log_lik[114]    -1.2e+00  3.2e-03    0.18   -1.5  -1.2e+00  -0.93     2915     3366      1.0
log_lik[115]    -1.5e+00  4.6e-03    0.25   -1.9  -1.4e+00   -1.1     2920     3371     1.00
log_lik[116]    -9.1e-01  1.6e-03   0.080   -1.1  -9.1e-01  -0.79     2589     2989     1.00

Samples were drawn using hmc with nuts.
For each parameter, N_Eff is a crude measure of effective sample size,
and R_hat is the potential scale reduction factor on split chains (at 
convergence, R_hat=1).
fit$summary()
# A tibble: 121 × 10
   varia…¹     mean   median     sd    mad      q5     q95  rhat ess_b…² ess_t…³
   <chr>      <dbl>    <dbl>  <dbl>  <dbl>   <dbl>   <dbl> <dbl>   <dbl>   <dbl>
 1 lp__    -5.33e+1 -5.29e+1 1.39   1.20   -55.9   -51.7    1.00   1678.   2278.
 2 beta[1] -3.89e-4 -1.80e-3 0.0882 0.0874  -0.148   0.147  1.00   3159.   2597.
 3 beta[2]  6.37e-1  6.42e-1 0.203  0.202    0.311   0.954  1.00   2407.   2351.
 4 beta[3] -7.36e-1 -7.38e-1 0.204  0.202   -1.07   -0.405  1.00   2471.   2540.
 5 sigma    9.64e-1  9.61e-1 0.0643 0.0640   0.863   1.07   1.00   2829.   2349.
 6 log_li… -1.05e+0 -1.03e+0 0.123  0.116   -1.27   -0.874  1.00   3739.   3523.
 7 log_li… -1.85e+0 -1.83e+0 0.286  0.286   -2.35   -1.43   1.00   4015.   3409.
 8 log_li… -9.35e-1 -9.30e-1 0.0815 0.0770  -1.08   -0.810  1.00   3014.   3257.
 9 log_li… -1.33e+0 -1.31e+0 0.167  0.163   -1.62   -1.09   1.00   3548.   3171.
10 log_li… -1.27e+0 -1.26e+0 0.149  0.145   -1.53   -1.06   1.00   3973.   2982.
# … with 111 more rows, and abbreviated variable names ¹​variable, ²​ess_bulk,
#   ³​ess_tail
# Extract draws for plotting.
draws <- fit$draws()
# Plot using bayesplot.
mcmc_hist(draws, pars = c("beta[1]", "beta[2]", "beta[3]"))
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
# Compute PSIS mode fit (p_loo).
loo(draws)
Warning: Relative effective sample sizes ('r_eff' argument) not specified.
For models fit with MCMC, the reported PSIS effective sample sizes and 
MCSE estimates will be over-optimistic.
Warning: Some Pareto k diagnostic values are too high. See help('pareto-k-diagnostic') for details.

Computed from 4000 by 121 log-likelihood matrix

         Estimate    SE
elpd_loo   -216.5  54.9
p_loo         6.8   3.0
looic       433.0 109.9
------
Monte Carlo SE of elpd_loo is NA.

Pareto k diagnostic values:
                         Count Pct.    Min. n_eff
(-Inf, 0.5]   (good)     120   99.2%   2293      
 (0.5, 0.7]   (ok)         0    0.0%   <NA>      
   (0.7, 1]   (bad)        0    0.0%   <NA>      
   (1, Inf)   (very bad)   1    0.8%   8         
See help('pareto-k-diagnostic') for details.
# Save the model output.
fit$save_object(file = here::here("Output", "foxes_fit.rds"))

Now when you want to load the model output into your report later, it’s simple.

fit <- readRDS(here::here("Output", "reedfrogs_fit.rds"))  # Base R.
fit <- read_rds(here::here("Output", "reedfrogs_fit.rds")) # Tidyverse.

Removing the Scaffolding

Stan has a number of awesome resources (including Statistical Rethinking): the Stan User’s Guide (plus the Stan Reference Manual and the Stan Functions Reference).

There are a growing number of supporting packages: {bayesplot}, {posterior}, {loo}, {ggdist}, {tidybayes}, and {brms}.